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SUMMARY 


This paper presents recent developments in advanced analysis methods for the computation of stress 
intensity factors and the J-integral under conditions of small scale yielding in structural panels with multi- 
site damage. The method of solution is based on the p-version of the finite element method. Its implemen- 
tation was designed to permit extraction of linear stress intensity factors using a superconvergent extraction 
method (known as the contour integral method) and evaluation of the J-integral following an elastic-plastic 
analysis. Coarse meshes are adequate for obtaining accurate results supported by p-convergence data. The 
elastic-plastic analysis is based on the deformation theory of plasticity and the von Mises yield criterion. 

The model problem consists of an aluminum plate with six equally spaced holes and a crack emanating 
from each hole. The cracks are of different sizes. The panel is subjected to a remote tensile load. Experi- 
mental results are available for the panel. The plasticity analysis provided the same limit load as the exper- 
imentally determined load. The results of elastic-plastic analysis were compared with the results of linear 
elastic analysis in an effort to evaluate how plastic zone sizes influence the crack growth rates. The onset of 
net-section yielding was determined also. The results show that crack growth rate is accelerated by the 
presence of adjacent damage, and the critical crack size is shorter when the effects of plasticity are taken 
into consideration. This work also addresses the effects of alternative stress-strain laws: The elastic-ideaily- 
plastic material model is compared against the Ramberg-Osgood model. 


INTRODUCTION 


Reliable stress intensity factor computation of cracked structural components is of major importance in 
modern design of aircraft structures where requirements for residual strength and fatigue crack propagation 
must be met In this paper the problem of computing the stress intensity factors for a structural panel with 
multi-site damage, with guarantee of reliability, is discussed. This work focuses on the accuracy and reli- 
ability of the numerical solution of a proposed mathematical model and on the influence that different mod- 
eling assumptions may have on the results of the analysis. Two kinds of error are important in this case: 

1. The differences between the exact solution of the mathematical problem formulated to represent a 
physical system or process and its numerical approximation are called errors of discretization. Is it 
possible to guarantee that the error of discretization is small? 
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. e differences between the exact solution of the mathematical problem formulated to represent a 
physical system and the actual response or behavior of the physical system are called errors of ideali- 
zation or modeling errors. How does incoiporating plasticity in the analysis affect the stress intensity 


In some cases the two errors, the errors of discretization and modeling, may partially cancel one another 
11 1S important to verify by means other than the experiment itself, that the numerical solution is 
close to the exact solution of the model. Only then is it possible to investigate whether the errors of ideali- 
zation are large or small by making comparisons with experimental observations. The efficient and reliable 
control of numerical errors, achieved by the use of a superconvergent method for the computation of the 

stress intensity factors, makes it feasible to investigate the sensitivity of crack extension to alternative mod- 
eling decisions. 

r ” * 1S P a P er > the eff ect of plasticity on the values of the fracture mechanics parameters is investigated, 
birst, the finite element implementation is discussed; and second the method is illustrated with the integrity 
assessment of a row of fastener holes with multi-site damage, for which experimental results are available 


FINITE ELEMENT IMPLEMENTATION 


Choice of an Extension Process 


In the finite element method the control of the errors of discretization can be achieved by mesh refine- 

r e ^ t / l enS1 ° n) : by mcreasin S the polynomial degree of elements (p-extension), or a combination of 
Doth (hp-extension). 

The size of a finite element is the diameter of the smallest circle (or sphere) that contains the element, 
is diameter is denoted by h and the diameter of the largest finite element in the mesh is denoted by h 
h-Extension involves letting -> 0. Alternatively, we can hold the number of elements constant and™** 
increase the polynomial degree of elements. The polynomial degree of elements is a vector p. p-Extension 
involves letting the smallest polynomial degree />„„■„-> i n hp-extensions mesh refinement is combined 
with an increase of the polynomial degree of elements. In these processes the number of equations that has 
to be solved, the number of degrees of freedom, is progressively increased, hence the name ‘extension’. 

ote that h- and p-extensions can be viewed as special cases of hp-extension, which is the general discreti- 
zation strategy of the finite element method. 

An important question is: “Which is the most efficient method of extension with respect to reducing dis- 
cretization errors . This question can be answered on the basis of a simple classification of the exact solu- 

categories- Pr ° blem 0 " e W ' SheS t0 Solve ‘ 1116 exact solutions ("ex) have been classified into three main 

Category A : u EX is analytic. A function is analytic at a point if it can be expanded into a Taylor series 
about that point on the entire solution domain, including the boundaries of the solution domain. Alterna- 
tively the domain can be divided into subdomains and u EX is analytic on each subdomain, including the 
boundary of each subdomain. The finite element mesh is so constructed that the boundaries of the subdo- 
mains are coincident with element boundaries. 

Category B: u EX is analytic on the entire domain, including the boundaries, with the exception of a finite 
number of points (in three dimensions u EX needs not be analytic along a finite number of lines) The 
mesh is so constructed that the points where u EX is not analytic are nodal points (in three dimensions the 



lines where u EX is not analytic are coincident with element edges). An example of problems in Cate- 
gory B is a cracked elastic body. The crack tip is a singular point 

Category C: Problems for which the exact solution is neither in Category A nor in Category B are in 
Category C. 

Most problems in mechanical and structural design belong in either Category A or Category B. Most 
problems in fluid dynamics belong in Category C. Knowing the classification of a problem, it is possible 
to compare h-, p- and hp-extensions on the basis of rate of convergence, that is, the rate of change of the 
error measured m energy norm with respect to the number of degrees of freedom. For problems in Cate- 
gory A and B the discretization errors are most efficiently controlled by p-extension. Designing a mesh 
such that when the polynomial degree of elements is uniformly increased then the discretization error 
decreases exponentially until the required degree of accuracy has been reached ; usually it involves grading 
the mesh in geometric progression in the vicinity of singular points (ref. 1). 


The Trunk and Product Spaces 


By 

rnl*i 


definition, the standard quadrilateral element Q® is the set of points (^,T|) which satisfy |£| <; 1, 
, see Figure 1. The trunk space of degree p is defined as the set of polynomials which can be 
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Figure 1. The standard quadrilateral element 


expressed as linear combinations of Q rj-i, i,j =0,1, ...p; i +j<p and supplemented by the monomial term £ 
T[ for p=l, and the monomial terms T|, £ yf for p > 2 . For example, for p=4 the trunk space includes all 
polynomials up to degree 4, plus two fifth order polynomials. The product space of degree p is the set of 
polynomials which can be expressed as linear combinations of Tj J , i,j =0,1, ...p. For example, for p=4 the 
product space includes all polynomials up to degree 8. 

The answer to the question: “Which space should be used in a given application?” depends on the 
accuracy requirements. The following points should be taken into consideration: 

1* The number of degrees of freedom increases much more rapidly with respect to p when the product 
space is used. For example, in the case of p=8 there are 47 degrees of freedom per field for the trunk 
space, and 8 1 degrees of freedom per field for the product space. 

2. For a given polynomial degree the error, measured in energy norm, is always smaller for the product 
space than for the trunk space. This is because the trunk space is a subset of the product space. 

3. For a given number of degrees of freedom the error for the product space, measured in energy norm, 
may be smaller or larger than the error for the trunk space. There appears to be no way of predicting 
which space performs better for specific cases. 
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. or static analyses, work performed by the computer increases more slowly with respect to increas- 
ing p for the product space than for the trunk space. This is because in the case of the product space 
e increase affects the internal modes only and those modes are eliminated by the solver locally i e 
no assembly is involved and the front width is unaffected. 

the I noT™nZl?Hr tati0n t> ^ hthe *""? and product spaces 216 available for quadrilateral elements, and 
the polynomial degree can be increased up to p=8. 


Elastic-Plastic Analysis 


an ma£erial n ° n ! inearities within ** Rework of the p-version was considered 

^ theasscss j ment of modeling assumptions. The elastic-plastic analysis is based on 
yon Mises yield criterion and the deformation theory of plasticity. The purpose and scope of this 
implementation are as follows (for further details see ref. 2). P 

1. Realistic mathematical models of real physical systems must have a capability to provide initial esti- 
mates for the effects of nonhneanties at a low computational cost. The deformation theory of plas- 
ticity serves this purpose well for a large class of practical problems. 

2. The effects of a single overload event on structures made of ductile materials are of substantial prac- 

ca importance. Such effects can be well represented by mathematical models based on the defor- 
mation theory of plasticity, provided that the plastic flow is contained, i.e., the plastic zone is 
surrounded by elastic material. 

3 ' P r °P a fi at jon of cracks in strain-hardening materials is generally correlated with the J-integral. 

ihe J-integral is based on the deformation theory of plasticity. 

4 - A n important feature of the linear implementation is that engineering data can be conveniently 
extracted from finite element solutions in the post-solution phase. The deformation theory of plastic- 
ity makes it feasible to extend this into the elastic-plastic regime because the data storage require- 
menis arc small. 

5 ' Jbf p ' version , 1S not susceptible to Poisson ratio locking and hence correct limit loads are obtained 
In the conventional (h-version) locking occurs when the displacement formulation is used For this 
reason alternative formulations, generally known as mixed methods, must be employed 


Elastic-Plastic Material Properties 


An important modeling assumption is the type of stress-strain law used for the elastic-plastic analysis. 
Four types of stress- strain relationships have been implemented: 

PoTson?S^f : t^ rial ^ haraaer ‘ 2e f by f° Ur Pieters (Figure 2): The modulus of elasticity (£), 
Poisson s ratio (tij the stress (S 70E ), which is the stress corresponding to the intersection of the stress 

^expmTiom 2 106 “ P8SSeS thr ° Ugh thC 0rigi " haS the slope of a 70 £ and an Ifponent „ in 


+ 


70 E 


(c— ) 


70 E 


( 1 ) 
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Typical values of n range between 4 and 90. Fractional values are permitted for n. In the limit as n — > oo, an 
elastic-perfectly plastic material behavior is approached. 

Elastoplastic: Elastic -perfectly plastic material characterized by E, v and the yield stress Sy 

Bilinear: Linear-elastic, linear strain hardening material characterized by E , t>, S y and the tangent modulus 
E t which characterizes strain hardening. y 

5-Parameter: Material characterized by two linear segments joined by a cubic spline. 

Solution Algorithm 


The following is an outline of the solution procedure implemented for two-dimensional elastic-plastic 
problems based on the deformation theory of plasticity. The procedure is known as direct integration. The 
iteration number is represented by a superscript in brackets: 

1. Obtain a linear solution for the problem. Ensure that the relative error in energy norm is small, cer- 
tainly under 5%, preferably under 1%. 


2. Compute the total equivalent strain £ in each Gauss point. The elastic equivalent strain is defined as: 

' 


l e - 




( 2 ) 


and the plastic equivalent strain is given as: 


IP = 



where x> is Poisson s ratio, £j, 82, £3 are the principal strains. The total equivalent strain is: 


(3) 


e = e e 4- eP . (4) 

3. Using (£)^ \ compute the secant modulus Ej k ^ at each Gauss point from the one-dimensional stress- 
strain curve. 


4. In each Gauss point for which £ > £ Y (the uniaxial strain at the onset of yielding), determine the elastic 
plastic material stiffness matrix using Ej^K Recompute the stiffness matrices for those elements for 
which £ > £ Y in one or more Gauss points, and obtain a new finite element solution u F J k+I K 
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5. Using and u FE f k+I K compute the strain and stress tensor components and the equivalent stress at 
each Gauss point. Using expressions (2)-(4), compute the total equivalent strain at each Gauss point 
( £ )(k+l) equivalent stress in each Gauss point satisfies: 


a+D 


<T|^ 


(5) 


then stop, else using (e)^ k+1 \ compute E s ^ +l \ increment kto k+1 and return to step 3. In (5), G r is the 
reference equivalent stress from the uniaxial stress- strain curve corresponding to (e)^ +1 \ and T is the 
pre-specified tolerance. The equivalent stress is computed based on the principal stresses as follows: 


o 



( 6 ) 


Fracture Mechanics Parameters 


In the neighborhood of the crack tip in an elastic body, the solution becomes singular and the stress val- 
ues go to infinity at the crack tip. Whether or not a crack will propagate, and at what rate, depends on the 
energy available to drive crack extension. In linear-elastic fracture mechanics the stress intensity factors are 
a measure of the stress singularity at the crack tip. They depend on the geometry of the body, the configura- 
tion of the crack, the boundary conditions and the loading. In the elastic-plastic regime, the J-integral has 
been increasingly used to characterize crack initiation and crack growth. J can be interpreted as the energy 
release rate, and in the linear-elastic case J=G (the energy release rate). 


Stress Intensity Factors 

The computation of Mode I and II stress intensity factors for linear elastic fracture mechanics was imple- 
mented using the Contour Integral Method, as described in ref. 1 : 

Kj = JtkA f V K u = J2ka[ 2) (7) 

where A ^ m) , m = 1, 2, is the first term in the asymptotic expansion of the solution in the neighborhood of 
the crack tip for modes I and II respectively. 

Let Fp be a circle of radius p centered on the crack tip (Figure 3), and assume that p is sufficiently close 
to the crack tip. Then we can write: 



where is an extraction function, T EE is the traction vector along Fp computed from the finite element 
solution u FE , and T w is the traction vector along Fp due to the extraction function. 


The J-integral 

The J-integral provides a mean to determine an energy release rate for cases where plasticity effects are 
not negligible. In the case of plane stress and plane strain the J-integral is defined by: 
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( 9 ) 


J = J {Wdy-T^ds) 

T x 

P 

where Fp is a contour around the crack tip, W is the strain energy density: 

e 

W - W (x, y) - ja.jdt.j. (10) 

0 

T is the traction vector (with components T x , T y ) along Fp, u is the displacement vector ( u x , u y ) and ds is 
an element of Fp. The coordinate system is located such that the origin is in the crack tip, the x-direction is 
parallel to the crack face (see Figure 3). 



If the integration is performed along a circular path of radius p centered at the crack tip, equation (9) can 
be written as: 

r \ a '\. 

J = J (™ 9 -^ ~ T yTx )(>dQ - di) 

—K 

The J-integral can be computed for linear and nonlinear isotropic or orthotropic materials. For the case 
of elastic-plastic materials, the integration path should be selected in such a way that it does not cut through 
the plastic zone around the crack tip. In the plastic region the strain energy density (W) is not defined. 


EXAMPLE PROBLEM 


Problem Description 


The model problem consists of a 2024-T3 aluminum alloy plate with six equally spaced holes and a 
crack emanating from each hole. The cracks are of different sizes and the panel is subjected to a remote ten- 
sile load. The panel was selected from ref. 3, and corresponds to an MSD test coupon for constant ampli- 
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tilde crack growth testing. Numerical and experimental crack growth data are available for the panel. 
Figure 4 shows the configuration of the panel for the example problem. 

The objective of the analysis is to compute the stress intensity factors and the J-integral for each crack at 
different stages during crack propagation. The results are compared with the experimental values. The elas- 
tic-plastic analysis is performed using two different stress-strain laws: Ramberg-Osgood and elastic-ide- 
ally-plastic. The length of each crack (measured from the edge of the hole) for the four cases analyzed is 


Material 2024-T3: 
Sy= 330 MPa 
E= 71 GPa 
0.33 

Dimensions (mm): 
e=12.5 
s=15.0 
D=5.0 

t=2.5 (thickness) 
W=100.0 
L=300.0 

Loading: 

a = 96.5 MPa 


Figure 4. Panel configuration for example problem. 

shown in Table 1. Case 1 corresponds to the initial crack length values; Case 2 represents the crack lengths 
after 6000 loading cycles; Case 3 after 12000 cycles; and Case 4 represents the crack lengths at failure 
(-13500 cycles). 


Table 1 : Crack Length (mm) for Each Case 
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Discretization 


In designing the finite element mesh and selecting the other discretization parameters, the smoothness of 
the exact solution of the mathematical model representing the physical system should be considered. For 
the example problem, the following modeling aspects are relevant: 

1. Two-dimensional (plane stress) analysis is considered appropriate for this case. 

2. Due to symmetry conditions, only half of the panel needs to be analyzed. 

3. Elastic-plastic material properties are used to assess the influence of crack tip plasticity on the fracture 
mechanics parameters. 

The solution of the mathematical model just described is smooth inside the domain and on the bound- 
aries, except at the crack tips where the stresses are unbounded in the linear case. For the nonlinear case 
when the material is elastic-perfectly plastic, the stresses are finite but the strains are singular at the crack 
tips. The design of the finite element mesh should account for the nature of the exact solution. Geometri- 
cally graded meshes toward the singular points are known to be optimal in conjunction with p-extension for 
problems of this type. However, the design of geometrically graded meshes for domains with multiple 
cracks emanating from holes with sizes comparable to the hole diameters is very cumbersome and in many 
cases impractical. For that reason, the use of very simple meshes combined with high order elements (prod- 
uct space) is selected. This strategy has been tested and documented in refs. 3 and 4. 

The finite element mesh consisting of 38 quadrilateral elements is shown in Figure 5. Each crack has 
been defined parametrically, so that the same mesh could be used for the four cases shown in Table 1 by 
adjusting the crack length parameters. Note that no mesh refinement around the crack tips has been used. 
Due to symmetry, only half the panel was discretized. 

Two elastic-plastic stress-strain relations were investigated. The Ramberg-Osgood material character- 
ized by: E=l\ GPa, 0=0.33, S 70E = 320 MPa, n= 12, and the elastic-perfectly plastic material characterized 
by S y = 330 MPa. The linear solution was computed for the product space and polynomial degree ranging 
from 1 to 8. Four nonlinear solutions were computed for each case in order to obtain an assessment of the 
discretization error for the nonlinear analysis. The linear runs corresponding to p-levels 5, 6, 7 and 8 were 
used to start the nonlinear iterations. The number of iterations required to complete any given analysis was 
based on a specified tolerance of 1.0% (see Eq. 5). 


Results 


Discretization Errors 

For each linear analysis, the estimated relative error in energy norm was computed from the eight avail- 
able solutions. Table 2 shows a typical result which corresponds to Case 1. The reported rate of conver- 
gence is typical for the p-extension in the case of category B type problems without geometrically graded 
meshes towards the singular points. For cases 1, 2 and 3, the estimated relative error in energy norm corre- 
sponding to p-leve!=5 and above was less than 5%. 
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Figure 5. Finite element mesh. 


Table 2: Estimated Relative Error in Energy Norm. Linear Solution for Case 1 


p-Level 

DOF 

Potential 
Energy [N mm] 

Rate of 
Convergence 

Estimated 
% Error j 

1 

117 

-2490.4413 

0.00 

10.23 

2 

392 

-2504.4379 

0.31 

7.02 

3 

819 

-2510.8236 

0.50 

4.86 

4 

1398 

-2513.3046 

0.51 

3.71 

5 

2129 

-2514.6060 

0.56 

2.93 

6 

3012 

-2515.3562 

0.62 

2.37 

7 

4047 

-2515.7996 

0.64 

1.96 

8 

5234 

-2516.0703 

0.64 

1.66 


The fracture mechanics parameters were computed for the available sequence of linear and nonlinear 
solutions: Eight runs for each linear case and 4 runs for each nonlinear case. Convergence of the extracted 
results was assessed in each case. Typical convergence data are shown in Table 3 for the linear and nonlin- 
ear fracture mechanics parameters of Case 2, Crack 1. 


Fracture Mechanics Parameters 

Tables 4 to 6 show the values of the fracture mechanics parameters (FMP) for each crack and for cases 1 
to 3 respectively. They include the stress intensity factor (Kj), the J-integral (J e ) computed from the linear 
solution and the values of the J-integral computed from the elastic-plastic solution for the Ramberg -Osgood 
material (J p R0 ) and the elastic-perfectly plastic material (J p EF> ). The relative difference between the elastic 





and plastic (Ramberg-Osgood) values of J is also included. All values included in the tables correspond to 
the linear or nonlinear solutions for p-level=8. 


Table 3: Convergence Data for Case 2, Crack 1 


p-Level 

K] [MPa m 1/2 ] 

J e [N rrf 1 ] 

J p RO [N m' 1 ] 

J p EP [N m-'j 

1 

11.26 

1148 


— 

2 

14.43 

2567 

— 

— 

3 

13.84 

2647 

— 

— 

4 

13.41 

2531 

— 

— 

5 

13.90 

2681 

2737 

2683 

6 

14.28 

2822 

2881 

2825 

7 

14.22 

2827 

2925 

2859 

8 

14.15 

2815 

2927 

2864 


Table 4: FMP for Case 1 (N=0) 


Parameter 

Value of the Parameter for Crack # 

1 

— 

2 

3 

4 

5 

6 


12.65 ' 

12.40 

12.49 

12.49 

12.11 

11.95 

JetNm- 1 ] 

2242 

2151 

2186 

2185 

2059 

2000 


2276 

2200 

2236 

2231 

2097 

2031 


2256 

2176 

2219 

2218 

2084 

2015 

| (Jp-J 'AW 

1.8 

2.3 

2.3 

1.8 

1.9 

1.5 


Table 5: FMP for Case 2 (N=60Q0) 





























































Table 6: FMP for Case 3 (N=12000) 


Parameter 

Value of the Parameter for Crack # 

1 

2 

3 

4 

5 

6 


18.78 

19.26 

19.04 

18.21 

16.67 

15.39 

J e [N m' 1 ] 

4885 

5152 

5019 

4604 

3893 

3325 


6292 

7057 

6681 

5786 

4504 

3603 


5585 

6072 

5771 

5151 

4196 

3497 

(J p RO -J c )/J e [%] 

28.6 

37.1 

33.1 

25.9 

15.7 

8.1 


These results clearly indicate that the stress intensity factors are not substantially affected by crack tip 
plasticity until the number of cycles is sufficiently high (near the end of the life of the panel). This is in 
agreement with the experimental observation: Towards the end of the life, the actual fatigue crack propaga- 
tion rate is larger than that predicted by using the linear stress intensity factors. Tables 7 to 9 show a com- 
parison of the experimental crack propagation rate (ref. 3) with those computed from the expression: 


= C(AK) m 
dN v y 


( 12 ) 


for the linear solution and from the expression: 


da 

dN 



(13) 


for the nonlinear solution. C and m are material constants (C=4.15 x 10 -8 , m- 2.2 when A K is given in ksi, 
and daldN is given in inch/cycle). The values of J p for the Ramberg-Osgood (R-O) and for the elastic-per- 
fectly-plastic material (E-P) are considered in Eq. 13. 


Table 7: Fatigue Crack Propagation Rates for Case 1 (N=0) 


Propagation Rate 
given by 
Eq. 12 
Eq. 13(R-Q) 
Eq. 13 (E-P) 
Experimental 


1 

2^28 

2.30 

2.28 

2.10 


2 

2T8 

2.21 

2.19 

2.20 


da/dN [10 4 mm/cycle] for Crack # 


3 


4 


2.21 

2.25 

2.23 

2.16 


2.21 

2.25 

2.23 

2.00 


5 


6 


2.07 

2.10 

2.09 

1.75 


2.01 

2.03 

2.01 

1.80 











































































Table 8: Fatigue Crack Propagation Rates for Case 2 (N=6000) 


Propagation Rate 


da/dN [10 4 mm/cycle] for Crack # 

given by 

a l 

a 2 

a 3 

M 

a 5 

36 

Eq. 12 

r 2.92 

I 2.93 

r 2.92 

2.86 

2.68 

2.49 

Eq. 13(R-0) 

3.03 

3.01 

3.00 

2.94 

r 2.72 

2.54 

Eq.13 (E-P) 

2.96 

2.93 

2.93 

2.87 

2.69 

2.50 

1 Experimental 

2.91 

2.69 

2.92 

2.64 

2.52 

2.20 


Table 9: Fatigue Crack Propagation Rates for Case 3 (N= 12000) 


Propagation Rate 


da/dN [10 4 mm/cycle] for Crack # 

given by 

a l 

a 2 

a 3 

a 4 

a 5 

36 

Eq. 12 

5.43 

5.74 

' 5.60 

5.08 

4.18 

j 3.49 

Eq. 13 (R-O) 

7.03 

7.98 

7.51 

6.41 

4.87 

3.81 

Eq. 13 (E-P) 

6.17 

6.76 

6.40 

5.64 

4.50 

3.69 

Experimental 

7.12 

8.95 

7.56 

6.40 

5.30 

3.47 


The Plastic Zone Size 


Figure 6 shows the extent of the plastic zone at each crack tip for cases 1, 2 and 3 for the Ramberg- 
Osgood stress-strain law. The plastic zone is the set of all points for which the total equivalent strain is 
greater than or equal to the uniaxial yield strain e > e Y , where Zy=Sy/E. 




Case 3: N=12000 


Figure 6. Plastic zone sizes for Ramberg-Osgood material. 
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Limit Load 


The limit load was computed for the elastic-perfectly plastic material, by imposing a uniform normal 
displacement along the top edge of the panel, and computing the stress resultant along that edge. The crack 
configuration of Case 4 (failure condition) was used, and a nonlinear solution was obtained for each value 
of the imposed displacement The results are shown in Figure 7. The limit load for this problem can be esti- 
mated as: 



Figure 7: Load-displacement relation for case 4, elastic-perfectly plastic material. 



(loads and constraints). This involves almost exclusively engineering time which will depend on the level 
of training of the analyst, and on the characteristics of the software tool used. The preparation of input data 
for the problem shown in Figure 5 required approximately 20 minutes of wall clock time. Since the crack 
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1. STRESS CHECK is a trademark of Engineering Software Research and Development, Inc., 7750 Clayton Rd., Suite 
204, St. Louis, MO 63117. 



lengths were defined in parametric form, and since the same description can be used to perform linear and 
nonlinear analyses, only one mesh was used for all cases analyzed. 

The execution time depends on the analysis type. To obtain a sequence of linear solutions (p-level=l to 
8), 1.5 minutes of CPU time were required for each case. For the nonlinear solutions, the CPU time 
depends on the number of iterations necessary to achieve the desired tolerance level (see Eq. 5) which in 
turn depends on the amount of plasticity associated with the solution. To obtain a sequence of nonlinear 

solutions (p-level-5 to 8) for the crack configuration of case 1, 9.5 minutes of CPU were required- for case 
2, 10.4 minutes of CPU; for case 3, 18.0 minutes of CPU. requires, tor case 

The extraction of the results (the fracture mechanics parameters) involves mostly analyst’s time. To 
compute the six crack tip stress intensity factors and J-integral (both elastic and plastic) for each case ana- 
lyzed, about 30 minutes of wall clock time were required. 

In summary, a complete analysis of a multiple-site damage panel for one crack configuration, including 
preparation of input data, execution of a sequence of Unear and nonlinear solutions, and extraction of the 
data ot interest, can be performed in about 1.5 hours of wall clock time. 

The Importance of Controlling Discretization Errors 

Intrinsic control of discretization errors by finite element programs is veiy important for properly corre- 
lating computed information with experimentally observed data and for using computed data in engineer- 
ing decisions which cannot be measured directly, and/or exploring a range of loadings and constraint 
conditions, not covered by the experiment. 

To illustrate this point, let <& EXI> be some experimental information, for example a measured displace- 
ment or strain, <D MOO the same information predicted by a mathematical model, for example the equations 
ot elasticity together with the appropriate material properties, boundary conditions and loading data; * 
the same information predicted by a finite element approximation of the mathematical model. <J> is gen- 
erally not known, nevertheless the purpose of an experiment is to determine whether the mathematical 
model correctly describes the physical system being modeled, that is, whether <J> MOD is sufficiently close to 
^exe* Writing: 


0 — 0 = (t> 

EXP ^ MOD ^ EXP 


0 - 4 - <J> — (t> < <f) 

^FEA^^FEA ^ MOD ~ V EXP 


d> + 0 - d> 

fea fea ^ mod 


(14) 


it is seen that if the quality of a mathematical model is to be assured then it is necessary to ensure that 
I “ ®mod\ 1S not larger, and preferably much smaller, than |0 £X/> - 0 FE J , which is known from mea- 
surement and computation. The quality control procedures of any finite element program should make it 
possible to ascertain that - <& Mop | is smaller than - . 

Unfortunately, there are many examples in industrial practice where close correlation with experimental 
results is achieved through skillful manipulation of the discretization parameters. The sensitivity of the 
approximate solution to discretization parameters is a strong indication that the numerical error is large, or 
the model is improperly defined. In the absence of intrinsic procedures designed for controlling approxima- 
tion errors, it is possible to correlate computed data with experimental observations through near cancella- 
tion of two large errors: 


0 


EXP 


0 „„ = 


FEA 


0 — 0 

EXP MOD 


±M 


+ 


0 — q> 

MOD ^ FEA 
— — — 

±M 


-0 


(15) 
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where M is some big number. Such models are unreliable and can be very misleading when the model is 
used for drawing general conclusions from experimental information. 


CONCLUSIONS 


Control of errors of discretization and modeling can be achieved if the analysis tools provide sufficient 
flexibility and are easy to use. It is also required that complex problems, both linear and nonlinear, be 
solved in hours instead of days. 

Ji . Advanced analysis methods for the computation of stress intensity factors and the J-integral under con- 
ditions of small scale plasticity have been implemented within the framework of the p-version of the finite 
element method. The implementation makes it possible for practicing engineers to assess damage in struc- 
tural components under elastic-plastic conditions with guaranteed reliability. The key elements which make 
this implementation unique can be summarized as follows: 

1. The use of the product space and superconvergent extraction for the stress intensity factors allows the 

use of very simple (coarse) meshes to obtain extremely reliable results for the linear elastic fracture 
mechanics parameters. 

2. The use of the deformation theory of plasticity with the von Mises yield criterion and the high-order 
finite elements (product space) makes it possible to use the same coarse meshes for the elastic-plastic 
analysis. Plastic zone sizes, J-integrals and limit loads can be easily and reliably computed from the 
finite element solutions. 

3. The availability of various representations for nonlinear material makes it possible to assess the influ- 
ence of different stress-strain laws in the results, 

4. The possibility of using spring boundary conditions, bearing loads, and precise geometric representa- 
tions allows for closer representation of actual problems. 

The application of these techniques to a multi-site damage panel clearly demonstrates that the investigation 
of different modeling assumptions can be accomplished in a reasonable amount of time, and that the results 
obtained can be used to assess the influence of plasticity on the residual strength of the panel. 
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